home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / debye.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  9.6 KB  |  388 lines
  1. /* specfunc/debye.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_debye.h>
  26.  
  27. #include "error.h"
  28. #include "check.h"
  29.  
  30. #include "chebyshev.h"
  31. #include "cheb_eval.c"
  32.  
  33. static double adeb1_data[17] = {
  34.    2.4006597190381410194,
  35.    0.1937213042189360089,
  36.   -0.62329124554895770e-02,
  37.    0.3511174770206480e-03,
  38.   -0.228222466701231e-04,
  39.    0.15805467875030e-05,
  40.   -0.1135378197072e-06,
  41.    0.83583361188e-08,
  42.   -0.6264424787e-09,
  43.    0.476033489e-10,
  44.   -0.36574154e-11,
  45.    0.2835431e-12,
  46.   -0.221473e-13,
  47.    0.17409e-14,
  48.   -0.1376e-15,
  49.    0.109e-16,
  50.   -0.9e-18
  51. };
  52. static cheb_series adeb1_cs = {
  53.   adeb1_data,
  54.   16,
  55.   -1.0, 1.0,
  56.   9
  57. };
  58.  
  59. static double adeb2_data[18] = {
  60.    2.5943810232570770282,
  61.    0.2863357204530719834,
  62.   -0.102062656158046713e-01,
  63.    0.6049109775346844e-03,
  64.   -0.405257658950210e-04,
  65.    0.28633826328811e-05,
  66.   -0.2086394303065e-06,
  67.    0.155237875826e-07,
  68.   -0.11731280087e-08,
  69.    0.897358589e-10,
  70.   -0.69317614e-11,
  71.    0.5398057e-12,
  72.   -0.423241e-13,
  73.    0.33378e-14,
  74.   -0.2645e-15,
  75.    0.211e-16,
  76.   -0.17e-17,
  77.    0.1e-18
  78. };
  79. static cheb_series adeb2_cs = {
  80.   adeb2_data,
  81.   17,
  82.   -1.0, 1.0,
  83.   10
  84. };
  85.  
  86. static double adeb3_data[17] = {
  87.    2.707737068327440945,
  88.    0.340068135211091751,
  89.   -0.12945150184440869e-01,
  90.    0.7963755380173816e-03,
  91.   -0.546360009590824e-04,
  92.    0.39243019598805e-05,
  93.   -0.2894032823539e-06,
  94.    0.217317613962e-07,
  95.   -0.16542099950e-08,
  96.    0.1272796189e-09,
  97.   -0.987963460e-11,
  98.    0.7725074e-12,
  99.   -0.607797e-13,
  100.    0.48076e-14,
  101.   -0.3820e-15,
  102.    0.305e-16,
  103.   -0.24e-17
  104. };
  105. static cheb_series adeb3_cs = {
  106.   adeb3_data,
  107.   16,
  108.   -1.0, 1.0,
  109.   10
  110. };
  111.  
  112. static double adeb4_data[17] = {
  113.    2.781869415020523460,
  114.    0.374976783526892863,
  115.   -0.14940907399031583e-01,
  116.    0.945679811437042e-03,
  117.   -0.66132916138933e-04,
  118.    0.4815632982144e-05,
  119.   -0.3588083958759e-06,
  120.    0.271601187416e-07,
  121.   -0.20807099122e-08,
  122.    0.1609383869e-09,
  123.   -0.125470979e-10,
  124.    0.9847265e-12,
  125.   -0.777237e-13,
  126.    0.61648e-14,
  127.   -0.4911e-15,
  128.    0.393e-16,
  129.   -0.32e-17
  130. };
  131. static cheb_series adeb4_cs = {
  132.   adeb4_data,
  133.   16,
  134.   -1.0, 1.0,
  135.   10
  136. };
  137.  
  138.  
  139. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  140.  
  141. int gsl_sf_debye_1_e(const double x, gsl_sf_result * result)
  142. {
  143.   const double val_infinity = 1.64493406684822644;
  144.   const double xcut = -GSL_LOG_DBL_MIN;
  145.  
  146.   /* CHECK_POINTER(result) */
  147.  
  148.   if(x < 0.0) {
  149.     DOMAIN_ERROR(result);
  150.   }
  151.   else if(x < 2.0*GSL_SQRT_DBL_EPSILON) {
  152.     result->val = 1.0 - 0.25*x + x*x/36.0;
  153.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  154.     return GSL_SUCCESS;
  155.   }
  156.   else if(x <= 4.0) {
  157.     const double t = x*x/8.0 - 1.0;
  158.     gsl_sf_result c;
  159.     cheb_eval_e(&adeb1_cs, t, &c);
  160.     result->val = c.val - 0.25 * x;
  161.     result->err = c.err + 0.25 * x * GSL_DBL_EPSILON;
  162.     return GSL_SUCCESS;
  163.   }
  164.   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
  165.     const int nexp = floor(xcut/x);
  166.     const double ex  = exp(-x);
  167.     double sum = 0.0;
  168.     double xk  = nexp * x;
  169.     double rk  = nexp;
  170.     int i;
  171.     for(i=nexp; i>=1; i--) {
  172.       sum *= ex;
  173.       sum += (1.0 + 1.0/xk)/rk;
  174.       rk -= 1.0;
  175.       xk -= x;
  176.     }
  177.     result->val = val_infinity/x - sum*ex;
  178.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  179.     return GSL_SUCCESS;
  180.   }
  181.   else if(x < xcut) {
  182.     result->val = (val_infinity - exp(-x)*(x+1.0)) / x;
  183.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  184.     return GSL_SUCCESS;
  185.   }
  186.   else {
  187.     result->val = val_infinity/x;
  188.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  189.     return GSL_SUCCESS;
  190.   }
  191. }
  192.  
  193.     
  194. int gsl_sf_debye_2_e(const double x, gsl_sf_result * result)
  195. {
  196.   const double val_infinity = 4.80822761263837714;
  197.   const double xcut = -GSL_LOG_DBL_MIN;
  198.  
  199.   /* CHECK_POINTER(result) */
  200.  
  201.   if(x < 0.0) {
  202.     DOMAIN_ERROR(result);
  203.   }
  204.   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
  205.     result->val = 1.0 - x/3.0 + x*x/24.0;
  206.     result->err = GSL_DBL_EPSILON * result->val;
  207.     return GSL_SUCCESS;
  208.   }
  209.   else if(x <= 4.0) {
  210.     const double t = x*x/8.0 - 1.0;
  211.     gsl_sf_result c;
  212.     cheb_eval_e(&adeb2_cs, t, &c);
  213.     result->val = c.val - x/3.0;
  214.     result->err = c.err + GSL_DBL_EPSILON * x/3.0;
  215.     return GSL_SUCCESS;
  216.   }
  217.   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
  218.     const int nexp = floor(xcut/x);
  219.     const double ex  = exp(-x);
  220.     double xk  = nexp * x;
  221.     double rk  = nexp;
  222.     double sum = 0.0;
  223.     int i;
  224.     for(i=nexp; i>=1; i--) {
  225.       sum *= ex;
  226.       sum += (1.0 + 2.0/xk + 2.0/(xk*xk)) / rk;
  227.       rk -= 1.0;
  228.       xk -= x;
  229.     }
  230.     result->val = val_infinity/(x*x) - 2.0 * sum * ex;
  231.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  232.     return GSL_SUCCESS;
  233.   }
  234.   else if(x < xcut) {
  235.     const double x2  = x*x;
  236.     const double sum = 2.0 + 2.0*x + x2;
  237.     result->val = (val_infinity - 2.0 * sum * exp(-x)) / x2;
  238.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  239.     return GSL_SUCCESS;
  240.   }
  241.   else {
  242.     result->val = (val_infinity/x)/x;
  243.     result->err = GSL_DBL_EPSILON * result->val;
  244.     CHECK_UNDERFLOW(result);
  245.     return GSL_SUCCESS;
  246.   }
  247. }
  248.  
  249.  
  250. int gsl_sf_debye_3_e(const double x, gsl_sf_result * result)
  251. {
  252.   const double val_infinity = 19.4818182068004875;
  253.   const double xcut = -GSL_LOG_DBL_MIN;
  254.  
  255.   /* CHECK_POINTER(result) */
  256.  
  257.   if(x < 0.0) {
  258.     DOMAIN_ERROR(result);
  259.   }
  260.   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
  261.     result->val = 1.0 - 3.0*x/8.0 + x*x/20.0;
  262.     result->err = GSL_DBL_EPSILON * result->val;
  263.     return GSL_SUCCESS;
  264.   }
  265.   else if(x <= 4.0) {
  266.     const double t = x*x/8.0 - 1.0;
  267.     gsl_sf_result c;
  268.     cheb_eval_e(&adeb3_cs, t, &c);
  269.     result->val = c.val - 0.375*x;
  270.     result->err = c.err + GSL_DBL_EPSILON * 0.375*x;
  271.     return GSL_SUCCESS;
  272.   }
  273.   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
  274.     const int nexp = floor(xcut/x);
  275.     const double ex  = exp(-x);
  276.     double xk  = nexp * x;
  277.     double rk  = nexp;
  278.     double sum = 0.0;
  279.     int i;
  280.     for(i=nexp; i>=1; i--) {
  281.       double xk_inv = 1.0/xk;
  282.       sum *= ex;
  283.       sum += (((6.0*xk_inv + 6.0)*xk_inv + 3.0)*xk_inv + 1.0) / rk;
  284.       rk -= 1.0;
  285.       xk -= x;
  286.     }
  287.     result->val = val_infinity/(x*x*x) - 3.0 * sum * ex;
  288.     result->err = GSL_DBL_EPSILON * result->val;
  289.     return GSL_SUCCESS;
  290.   }
  291.   else if(x < xcut) {
  292.     const double x3 = x*x*x;
  293.     const double sum = 6.0 + 6.0*x + 3.0*x*x + x3;
  294.     result->val = (val_infinity - 3.0 * sum * exp(-x)) / x3;
  295.     result->err = GSL_DBL_EPSILON * result->val;
  296.     return GSL_SUCCESS;
  297.   }
  298.   else {
  299.     result->val = ((val_infinity/x)/x)/x;
  300.     result->err = GSL_DBL_EPSILON * result->val;
  301.     CHECK_UNDERFLOW(result);
  302.     return GSL_SUCCESS;
  303.   }
  304. }
  305.  
  306.  
  307. int gsl_sf_debye_4_e(const double x, gsl_sf_result * result)
  308. {
  309.   const double val_infinity = 99.5450644937635129;
  310.   const double xcut = -GSL_LOG_DBL_MIN;
  311.  
  312.   /* CHECK_POINTER(result) */
  313.  
  314.   if(x < 0.0) {
  315.     DOMAIN_ERROR(result);
  316.   }
  317.   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
  318.     result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;
  319.     result->err = GSL_DBL_EPSILON * result->val;
  320.     return GSL_SUCCESS;
  321.   }
  322.   else if(x <= 4.0) {
  323.     const double t = x*x/8.0 - 1.0;
  324.     gsl_sf_result c;
  325.     cheb_eval_e(&adeb4_cs, t, &c);
  326.     result->val = c.val - 2.0*x/5.0;
  327.     result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;
  328.     return GSL_SUCCESS;
  329.   }
  330.   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
  331.     const int nexp = floor(xcut/x);
  332.     const double ex  = exp(-x);
  333.     double xk  = nexp * x;
  334.     double rk  = nexp;
  335.     double sum = 0.0;
  336.     int i;
  337.     for(i=nexp; i>=1; i--) {
  338.       double xk_inv = 1.0/xk;
  339.       sum *= ex;
  340.       sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;
  341.       rk -= 1.0;
  342.       xk -= x;
  343.     }
  344.     result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;
  345.     result->err = GSL_DBL_EPSILON * result->val;
  346.     return GSL_SUCCESS;
  347.   }
  348.   else if(x < xcut) {
  349.     const double x2 = x*x;
  350.     const double x4 = x2*x2;
  351.     const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;
  352.     result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;
  353.     result->err = GSL_DBL_EPSILON * result->val;
  354.     return GSL_SUCCESS;
  355.   }
  356.   else {
  357.     result->val = (((val_infinity/x)/x)/x)/x;
  358.     result->err = GSL_DBL_EPSILON * result->val;
  359.     CHECK_UNDERFLOW(result);
  360.     return GSL_SUCCESS;
  361.   }
  362. }
  363.  
  364.  
  365. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  366.  
  367. #include "eval.h"
  368.  
  369. double gsl_sf_debye_1(const double x)
  370. {
  371.   EVAL_RESULT(gsl_sf_debye_1_e(x, &result));
  372. }
  373.  
  374. double gsl_sf_debye_2(const double x)
  375. {
  376.   EVAL_RESULT(gsl_sf_debye_2_e(x, &result));
  377. }
  378.  
  379. double gsl_sf_debye_3(const double x)
  380. {
  381.   EVAL_RESULT(gsl_sf_debye_3_e(x, &result));
  382. }
  383.  
  384. double gsl_sf_debye_4(const double x)
  385. {
  386.   EVAL_RESULT(gsl_sf_debye_4_e(x, &result));
  387. }
  388.